##### Data format ##### -----
rm(list=ls())
setwd("D:/OneDrive - zju.edu.cn/data/public_database/cBioportal/prad_su2c_2019/prad_su2c_2019/")

data_raw <- read.table("data_mRNA_seq_fpkm_polya.txt",header=T)
data <- as.matrix(data_raw[,-1])
row.names(data) <- data_raw[,1]
tmp <- by(data,
          rownames(data),
          function(x)x[which.max(rowMeans(x)),])
tmp1 <- do.call(rbind,tmp)
data <- tmp1
data <- as.matrix(data)
write.table(data,"expression_FPKM.txt",sep="\t")

data = read.table("expression_FPKM.txt",header=T,sep="\t")

patient = read.table("data_clinical_patient.txt",header = T,sep = "\t")
sample = read.table("data_clinical_sample.txt",header = T,sep = "\t")
library(stringr)
library(dplyr)

meta = left_join(sample,patient,by="PATIENT_ID")
meta$SAMPLE_ID = gsub("\\_","\\.",meta$SAMPLE_ID,ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
meta$SAMPLE_ID = gsub("\\-","\\.",meta$SAMPLE_ID,ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
rownames(meta) = meta$SAMPLE_ID

colnames(data) = gsub("\\_","\\.",colnames(data),ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
colnames(data) = gsub("\\-","\\.",colnames(data),ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
meta_expr = meta[colnames(data),]

write.table(data,"expression_FPKM.txt",sep="\t")
write.table(meta_expr,"expression_meta.txt",sep="\t")

meta_UNK = meta_expr[meta_expr$ABI_ENZA_EXPOSURE_STATUS == "UNK",]
meta_OnT = meta_expr[meta_expr$ABI_ENZA_EXPOSURE_STATUS == "On treatment",]
meta_Naive = meta_expr[meta_expr$ABI_ENZA_EXPOSURE_STATUS == "Naive",]
meta_Exposed = meta_expr[meta_expr$ABI_ENZA_EXPOSURE_STATUS == "Exposed",]

gene = c("GREM1", "FGF1", "FGF2", "FGF3", "FGF4", "FGF5", "FGF6", "FGF7", "FGF8", "FGF9",
         "FGF10", "FGF16", "FGF17", "FGF18", "FGF19", "FGF20", "FGF21", "FGF22", "FGF23")

data = log2(data+1)

expr_UNK = data[gene,rownames(meta_UNK)]
expr_OnT = data[gene,rownames(meta_OnT)]
expr_Naive = data[gene,rownames(meta_Naive)]
expr_Exposed = data[gene,rownames(meta_Exposed)]

write.table(expr_UNK,"expr_UNK.txt",sep="\t")
write.table(expr_OnT,"expr_OnT.txt",sep="\t")
write.table(expr_Naive,"expr_Naive.txt",sep="\t")
write.table(expr_Exposed,"expr_Exposed.txt",sep="\t")
##### GREM1 & AR score ##### -----
### data format ###
setwd("D:/OneDrive - zju.edu.cn/data/public_database/cBioportal/prad_su2c_2019/prad_su2c_2019/")

exprset = read.table("data_mRNA_seq_fpkm_polya.txt",header=T,sep="\t")
matrix = as.matrix(exprset[,2:ncol(exprset)])
rownames(matrix) = exprset[,1]

meta = read.table("data_clinical_sample.txt",sep="\t",header=T)
meta$SAMPLE_ID = gsub("\\_","\\.",meta$SAMPLE_ID,ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
meta$SAMPLE_ID = gsub("\\-","\\.",meta$SAMPLE_ID,ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
rownames(meta) = meta$SAMPLE_ID

tmp = by(matrix,rownames(matrix),function(x)x[which.max(rowMeans(x)),])
exprset = do.call(rbind,tmp)
colnames(exprset) = gsub("\\_","\\.",colnames(exprset),ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
colnames(exprset) = gsub("\\-","\\.",colnames(exprset),ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
write.csv(exprset,"expression_FPKM_polya.txt")
exprset = read.csv("expression_FPKM_polya.txt",header=T,row.names = 1)

rna_meta = meta[colnames(exprset),]
write.csv(rna_meta,"rna_meta.txt")
gene = "GREM1"



### ARSi ###
meta_ARSI = read.table("ARSi_time.txt",sep="\t",header=T)
colnames(meta_ARSI)[1] = "PATIENT_ID"
rownames(meta_ARSI) = meta_ARSI$PATIENT_ID
meta_patient = read.table("data_clinical_patient.txt",header = T,sep = "\t")
rownames(meta_patient) = meta_patient$PATIENT_ID
meta_patient = meta_patient[rownames(meta_ARSI),]

meta = read.table("data_clinical_sample.txt",sep="\t",header=T)
meta$SAMPLE_ID = gsub("\\_","\\.",meta$SAMPLE_ID,ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
meta$SAMPLE_ID = gsub("\\-","\\.",meta$SAMPLE_ID,ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
rownames(meta) = meta$SAMPLE_ID
meta_sample = meta[rownames(FGF),]

meta_ARSI = merge(meta_ARSI,meta_patient,by="PATIENT_ID")
meta_ARSI = merge(meta_ARSI,meta_sample,by="PATIENT_ID")
Grem1_ARSI = Grem1_expr[meta_ARSI$Sample.ID,]
meta_Grem1_ARSI = cbind(Grem1_ARSI,meta_ARSI)
meta_Grem1_ARSI = na.omit(meta_Grem1_ARSI[,1:6])
write.csv(meta_Grem1_ARSI,"meta_Grem1_ARSI.csv")

write.csv(meta_sample,"meta_sample.csv")
write.csv(meta_ARSI,"meta_ARSI.csv")

median_Grem1 = median(meta_Grem1_ARSI$GREM1)

high_group = meta_Grem1_ARSI[meta_Grem1_ARSI$GREM1>median_Grem1,]
high_group$group = "high"
high_group = high_group[order(high_group$GREM1),]
#high_group = high_group[(((nrow(high_group)/2)+1):nrow(high_group)),]
low_group = meta_Grem1_ARSI[meta_Grem1_ARSI$GREM1<median_Grem1,]
low_group$group = "low"
low_group = low_group[order(low_group$GREM1),]
#low_group = low_group[1:(nrow(low_group)/2),]
dat = rbind(high_group,low_group)
library(survival)
library(survminer)

### survival curve ###
my.surv <- Surv(dat$OS.from.first.line.ARSI.start,dat$Death.status=='1')
kmfit1 <- survfit(my.surv~group,data = dat)
ggsurvplot(kmfit1,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE,legend.title="",xlab = "Months")

### c index ###
sum.surv = summary(coxph(Surv(dat$OS.from.first.line.ARSI.start,dat$Death.status=='1')~group,data=dat))
c_index_se = sum.surv$concordance
c_index = c_index_se[1]
c_index.ci_low = c_index - c_index_se[2]
c_index.ci_high = c_index + c_index_se[2]



##### GREM1 & FGFs ##### -----
setwd("D:/OneDrive - zju.edu.cn/data/public_database/cBioportal/prad_su2c_2019/prad_su2c_2019/")

exprset = read.table("GREM1_and_FGFs.txt",sep="\t",header=T,row.names = 1)

library(RColorBrewer)
library(pheatmap)
color <- rev(colorRampPalette(brewer.pal(10,"RdBu"))(10))

expr_zscore <- t(scale(t(exprset)))
expr_zscore[expr_zscore > 2] <- 2
expr_zscore[expr_zscore < -2] <- -2

pheatmap::pheatmap(expr_zscore, cluster_cols = T, cluster_rows = F)

data = read.table("expression_FPKM.txt",header=T,sep="\t")
library(GSVA)
list = data.frame(FGF_gene=c("FGF1","FGF2","FGF3","FGF4","FGF5","FGF6","FGF7","FGF8","FGF9",
                            "FGF10","FGF16","FGF17","FGF18","FGF19","FGF20","FGF21","FGF22","FGF23"))
data = log2(data+1)
FGF_score = as.data.frame(t(gsva(as.matrix(data),list2)))
GREM1 = t(data["GREM1",])

cor = cbind(GREM1,FGF_score)
write.csv(cor,"GREM1_AR_score.csv")

##### GREM1 & FGFs correlation ##### -----
setwd("D:/OneDrive - zju.edu.cn/data/public_database/cBioportal/prad_su2c_2019/prad_su2c_2019/")
#exprset = read.table("GREM1_and_FGFs.txt",sep="\t",header=T,row.names = 1)

data = read.table("expression_FPKM.txt",header=T,sep="\t")
FGF_gene=c("FGF1","FGF2","FGF3","FGF4","FGF5","FGF6","FGF7","FGF8","FGF9",
           "FGF10","FGF16","FGF17","FGF18","FGF19","FGF20","FGF21","FGF22","FGF23")
Grem1 = "GREM1"
data = log2(data+1)

### total ###
exprset = t(data[c(Grem1,FGF_gene),])
write.csv(exprset,"GREM1_FGF_expression.csv")

exprset = data[c(Grem1,FGF_gene),]
cor <- cor(t(exprset))

pdf("Grem1_FGFs_correlation.pdf")
corrplot::corrplot(-cor, method = "shade", type = "full", shade.col = NA, tl.col = "black",tl.cex=0.5)
dev.off()


##### AR_score ##### -----
data = read.table("expression_FPKM.txt",header=T,sep="\t")
data = log2(data+1)

library(GSVA)

list1 = data.frame(AR_gene=c("AR","KLK3","KLK2","TMPRSS2","NKX3-1","PLPP1","PMEPA1","ALDH1A3","STEAP4","FKBP5"))
list2 = data.frame(AR_gene=c("KLK3","KLK2","TMPRSS2","NKX3-1","GNMT","PMEPA1","MPHOSPH9","ZBTB10","EAF2","CENPN",
                             "C1orf116","ACSL3","PTGER4","ABCC4","NNMT","ADAM7","ELL2","MED28","HERC3","MAF","TNK1",
                             "GLRA2","MAPRE2","PIP4K2B","MAN1A1","CD200","FKBP5"))
list3 = data.frame(NE_gene=c("SYP","CHGA","CHGB","ENO2","CHRNB2","SCG3","SCN3A","PCSK1","ELAVL4","NKX2-1"))

### peter nelson ###
AR_score = as.data.frame(t(gsva(as.matrix(data),list1)))
ARPPC = data[,AR_score$AR_gene>0]
ARNPC = data[,AR_score$AR_gene<0]

#Grem1_ARPPC = log2(t(ARPPC["GREM1",])+1)
#Grem1_ARNPC = log2(t(ARNPC["GREM1",])+1)
Grem1_ARPPC = t(ARPPC["GREM1",])
Grem1_ARNPC = t(ARNPC["GREM1",])

write.csv(Grem1_ARPPC,"Grem1_ARPPC.csv")
write.csv(Grem1_ARNPC,"Grem1_ARNPC.csv")

### Hieronymus 2006 ###
AR_score = as.data.frame(t(gsva(as.matrix(data),list2)))
ARPPC = data[,AR_score$AR_gene>0]
ARNPC = data[,AR_score$AR_gene<0]

#Grem1_ARPPC = log2(t(ARPPC["GREM1",])+1)
#Grem1_ARNPC = log2(t(ARNPC["GREM1",])+1)
Grem1_ARPPC = t(ARPPC["GREM1",])
ARPPC_score = as.data.frame(AR_score[colnames(ARPPC),])
colnames(ARPPC_score) = "AR_score"
rownames(ARPPC_score) = colnames(ARPPC)
Grem1_ARPPC = cbind(Grem1_ARPPC,ARPPC_score)

Grem1_ARNPC = t(ARNPC["GREM1",])
ARNPC_score = as.data.frame(AR_score[colnames(ARNPC),])
colnames(ARNPC_score) = "AR_score"
rownames(ARNPC_score) = colnames(ARNPC)
Grem1_ARNPC = cbind(Grem1_ARNPC,ARNPC_score)

write.csv(Grem1_ARPPC,"Grem1_ARPPC.csv")
write.csv(Grem1_ARNPC,"Grem1_ARNPC.csv")

##### ARSi time ##### -----
ARSi = read.table("ARSi_time.txt",header=T,sep="\t")
colnames(ARSi)[1] = "PATIENT_ID"



##### correlation #####
library(GSVA)
library(corrplot)
library(psych)
library(Hmisc)
library(tidyr)
library(tidyverse)
library(ggplot2)
library(RColorBrewer)
color <- rev(colorRampPalette(brewer.pal(10,"RdBu"))(10))

exprset = read.csv("expression_FPKM_polya.txt",header=T,row.names = 1)
rna_meta = read.csv("rna_meta.txt",header=T,row.names = 1)

AR_score = as.data.frame(rna_meta$AR_SCORE)
AR_score_1 = AR_score
data = log2(exprset+1)
list2 = data.frame(AR_gene=c("AR","KLK3","KLK2","TMPRSS2","NKX3-1","PLPP1","PMEPA1","ALDH1A3","STEAP4","FKBP5"))
AR_score_2 = as.data.frame(t(gsva(as.matrix(data),list2)))
colnames(AR_score_1) = "AR_score_1"
colnames(AR_score_2) = "AR_score_2"

list2 = as.character(list2$AR_gene)
expr_cor = t(data[c("GREM1",list2),])
expr_cor_1 = as.matrix(na.omit(cbind(AR_score_1,AR_score_2,expr_cor)))

cor_2 = cor(expr_cor_1)
corrplot::corrplot(-cor_2, method = "shade", type = "full", shade.col = NA, tl.col = "black",tl.cex=0.5)


cor_1 = rcorr(expr_cor_1)

cor_1_r = as.data.frame(cor_1$r)
cor_1_p = as.data.frame(cor_1$P)

cor_1_p[is.na(cor_1_p)] <-1

cor_1_p=data.frame(gene1=rownames(cor_1_p),cor_1_p)
cor_1_r=data.frame(gene1=rownames(cor_1_r),cor_1_r)

data1 <- cor_1_p %>% 
  pivot_longer(cols=-1,
               names_to= "gene2",
               values_to = "p")
data2 <- cor_1_r %>% 
  pivot_longer(cols=-1,
               names_to= "gene2",
               values_to = "r")

data1$r=data2$r
data1=as.data.frame(data1)
data1$gene1 = as.character(data1$gene1)

data1[data1=="AR_score_1"]<-"A"
data1[data1=="AR_score_2"]<-"B"
data1[data1=="GREM1"]<-"C"
data1[data1=="AR"]<-"D"
data1[data1=="KLK3"]<-"E"
data1[data1=="KLK2"]<-"F"
data1[data1=="TMPRSS2"]<-"G"
data1[data1=="NKX3.1"]<-"H"
data1[data1=="NKX3-1"]<-"H"
data1[data1=="PLPP1"]<-"I"
data1[data1=="PMEPA1"]<-"J"
data1[data1=="ALDH1A3"]<-"K"
data1[data1=="STEAP4"]<-"L"
data1[data1=="FKBP5"]<-"M"

data1$p[data1$p>0.05]<-1
data1$p[data1$p<=0.05&data1$p>0.01]<-2
data1$p[data1$p<=0.01&data1$p>0.001]<-3
data1$p[data1$p<=0.001&data1$p>=0]<-4


p1<-ggplot(data1,aes(gene1,gene2))+
  geom_point(shape=21,aes(fill=r,size=p))+
  theme_minimal()+xlab(NULL)+ylab(NULL)+
  scale_size_continuous(range=c(7,15))+
  scale_fill_gradientn(colours=color,limits=c(-1,1)) +
  theme(legend.position = "left",legend.box = "vertical",panel.grid.major =element_blank(),panel.border = element_rect(fill=NA,color="black", size=0.6, linetype="solid"),axis.ticks = element_line(size = 0.6),
        axis.text.x = element_text(angle = 20))+
  theme(aspect.ratio=1)
p1





  
list3 = data.frame(AR_gene=c("AR","PLPP1","ALDH1A3","STEAP4","KLK3","KLK2","TMPRSS2","NKX3-1","GNMT","PMEPA1","MPHOSPH9","ZBTB10","EAF2","CENPN",
                               "C1orf116","ACSL3","PTGER4","ABCC4","NNMT","ADAM7","ELL2","MED28","HERC3","MAF","TNK1",
                               "GLRA2","MAPRE2","PIP4K2B","MAN1A1","CD200","FKBP5"))
list3 = as.character(list3$AR_gene)